(12) INTERNATIONAL APPLICATION PUBLISHED UNDER THE PATENT COOPERATION TREATY (PCT) 



(19) World Intellectual Property 
Organization 
International Bureau 

(43) Internationa) Publication Date 
7 April 2005 (07.04.2005) 




PCT 



III 1 1] 11 1 1 II II 1 1 I! II 1 1 II II II III! II I II 

r 

(10) International Publication Number 

WO 2005/031630 A2 



(51) International Patent Classification 7 : G06F 19/00 

(21) International Application Number: 

PCT/US 2004/03 0571 

(22) International Filing Date: 

17 September 2004 (17.09.2004) 



(25) Filing Language: 

(26) Publication Language: 



English 
English 



(30) Priority Data: 

10/667,045 22 September 2003 (22.09.2003) US 

(71) Applicant: UT-BATTELLE LLC [US/US]; Oak Ridge 
National Laboratory, Post Office Box 2008, Oak Ridge, TN 
37831 (US). 

(72) Inventor: HTVELY, Lee, M.; 192 Williams Road, 
Philadelphia, TN 37846 (US). 

(74) Agent: MCGOVERN, Michael, J.; Quarles & 
Brady LLP, 411 E. Wisconsin Avenue, Milwaukee, 
WI 53202^497 (US). 



(81) Designated States (unless otherwise indicated, for every 
kind of national protection available): AE, AG, AL, AM, 
AT, AU, AZ, BA, BB, BG, BR, BW, BY, BZ, CA, CH, CN, 
CO, CR, CU, CZ, DE, DK, DM, DZ, EC, EE, EG, ES, FI, 
GB, GD, GE, GH, GM, HR, HU, ID, IL, IN, IS, JP, KE, 
KG, KP, KR, KZ, LC, LK, LR, LS, LT, LU, LV, MA, MD, 
MG, MK, MN, MW, MX, MZ, NA, NI, NO, NZ, OM, PG, 
PH, PL, PT, RO, RU, SC, SD, SE, SG, SK, SL, SY. TJ, TM. 
TN, TR, TT, TZ, UA, UG, US, UZ, VC, VN, YU, ZA, ZM, 
ZW. 

(84) Designated States (unless otherwise indicated, for every 
kind of regional protection available): ARIPO (BW, GH, 
GM, KE, LS, MW, MZ, NA, SD, SL, SZ, TZ, UG, ZM, 
ZW), Eurasian (AM, AZ, BY, KG, KZ, MD, RU, TJ, TM), 
European (AT, BE, BG, CH, CY, CZ, DE, DK, EE, ES, FI, 
FR, GB, GR, HU, IE, IT, LU, MC, NL, PL, PT, RO, SE, SI, 
SK, TR), OAPI (BF, BJ, CF, CG, CI, CM, GA, GN, GQ, 
GW, ML, MR, NE, SN, TD, TG). 

Published: 

— without international search report and to be republished 
upon receipt of that report 

For two-letter codes and other abbreviations, refer to the "Guid- 
ance Notes on Codes and Abbreviations" appearing at the begin- 
ning of each regular issue of the PCT Gazette. 



(54) Title: METHODS FOR IMPROVED FOREWARNING OF CRITICAL EVENTS ACROSS MULTIPLE DATA CHANNELS 



^ 2 L< 8 > 
o 



17 



-J U 




O (57) Abstract: This disclosed invention concerns improvements in forewarning of critical events via phase-space dissimilarity anal- 
ly ysis of data from mechanical devices, electrical devices, biomedical data, and other physical processes. First, a single channel of 
O process-indicative data is selected that can be used in place of multiple data channels without sacrificing consistent forewarning of 
O critical events. Second, the method discards data of inadequate quality via statistical analysis of the raw data, because the analysis 

of poor quality data always yields inferior results. Third, two separate filtering operations are used in sequence to remove both 
Q high-frequency and low-frequency artifacts using a zero-phase quadratic filter. Fourth, the method constructs phase-space dissimi- 
^ larity measures (PSDM) by combining of multi-channel time-serial data into a multi-channel time-delay phase-space reconstruction. 

Fifth, the method uses a composite measure of dissimilarity (C,) to provide a forewarning of failure and an indicator of failure onset 
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METHODS FOR IMPROVED FOREWARNING OF CRITICAL EVENTS 
ACROSS MULTIPLE DATA CHANNELS 



BACKGROUND OF THE INVENTION 



[0001] The field of the invention is computer methods for 
analyzing nonlinear processes to forewarn of critical events 
in nonlinear processes. 

[0002] Examples of critical events are mechanical or 
electrical failures in machines; epileptic seizures, 
ventricular fibrillations, fainting (syncope), breathing 
difficulty, and sepsis in human medical patients; and other 
physical processes. Further examples of nonlinear processes 
include brain waves, heart waves, chest sounds, transients in 
power systems, airflow over automobiles and airplanes, weather 
and climate dynamics, water flow around submarines, machine 
tool -part interaction (e.g., tool chatter), nuclear reactor 
instabilities, fusion plasma instabilities, earthquakes, 
turbulent flow in pipes, planetary/satellite motion. 

[0003] Engineering, medical, and research applications 
frequently must distinguish a difference between two 
apparently similar but actually different states in a 
nonlinear process. Examples include various data for machine 
failure, pre-seizure versus non-seizure brain waves, pre- 
fibrillation versus fibrillation heart waves, pre-syncope 
versus syncope heart waves, pre-sepsis versus sepsis heart 
waves, and normal versus abnormal chest sounds as an indicator 
of breathing difficulty. The electrical/mechanical community 
calls this problem, "condition monitoring." The equivalent term 
in the medical community, is diagnostic, medical, or health 
monitoring. In the computer/networking world, the term is 
"security monitoring . " 

[0004] Hively et al., U.S. Pat. Nos. 5,743,860 and 5,857,978 
disclose methods for detecting and predicting epileptic 
seizures by acquiring brain wave data from a patient, and 
analyzing the data with traditional nonlinear methods. 
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[0005] Hively et al., U.S. Pat. No. 5,815,413 discloses the 
applicability of nonlinear techniques to monitor machine 
conditions, such as the condition of a drill bit or the 
performance of an electrical motor driving a pump. 

[0006] Clapp et al., U.S. Pat. No. 5,626,145 discloses the 
removal of artifacts representing eye blinks from EEG data 
using a zero-phase . quadratic filter. As one normally skilled 
in the art can appreciate, the same method can remove other 
artifacts from data, such as quasi -periodic variations from 
three-phase electrical current, voltage, or power; and quasi- 
periodic oscillations from one or more channels of 
acceleration; and breathing artifacts from heart waves and 
chest sounds ; . 

[0007] Hively et' al., U.S. Pat. No. 6,484,132 discloses a 
plurality of nonlinear techniques for using and enhancing 
..phase space dissimilarity measures (PSDM) to forewarn of 
machine failures, as well as impending epileptic events from 
scalp EEG in ambulatory settings. PSDM yield superior 
performance over traditional nonlinear indicators, such as 
Kolmogorov entropy, Lyapunov exponents, and correlation 
dimension. 

SUMMARY OF THE INVENTION 

[0008] The present invention improves on prior disclosed 
methods involving the use of phase space dissimilarity 
measures (PSDM) to provide forewarning indications of critical 
events. Such events include failures in critical mechanical 
or electrical equipment, in medical treatment of human 
patients, and in many other applications. 

[0009] The invention provides for better selection of 
process -indicative data for such analysis, such that one data 
channel can be used in place of multiple data channels without 
sacrificing consistent forewarning of critical events. 
Examples of such a process-indicative data are vibration power 
for mechanical systems, electrical power for electrically- 
driven machines, and differences between adjacent scalp EEG 
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(electroencephalogram) channels. 

[0010] The invention also provides a novel method to check 
the quality of the data, such that batches of inadequate 
quality data can be identified and replaced by adequate 
quality data. 

[0011] The invention also provides a novel method including 
at least two separate filtering operations to remove both 
high-frequency and low- frequency artifacts using a zero-phase 
quadratic filter, which previously has been used only removing 
low-frequency artifacts in the prior art. 

[0012] The invention also gives a novel method to combine 
several channels of time-serial data into a dC-dimensional 
multi-channel phase-space (PS) vector, which has the form: 

s i= [sa; i( s(i)^ , ~, B(u Ma .„x, ~, s(o it sfc;^- b (o 1+ftM ,d . 

Here, s £ {j) denotes the symbolized data for j-th channel at 
time t,. C denotes the total number of data channels in the 
nrulti -channel phase -space reconstruction. 

[0013] The invention also provides an end-of-life 
forewarning indicator (G) , which is computed from a composite 
measure of dissimilarity (C) . Here C, is the sum of the four 
renormalized PSDMs for the i-th subset: 



(1) 



The end-of-life indication (G) is quantified from this 
composite measure as follows. The appropriate process- 
indicative data is selected and acquired (for example) over 
10-second snapshots at periodic intervals (such as once per 
minute), and then divided into subsets of N data points, which 
typically has a value in the range, 10,000 < N Z 100,000. The 
analysis then proceeds as follows: (a) each subset is artifact 
filtered; (b) the data in each subset is symbolized; (c) each 
subset is converted into a sequence of multi-channel phase- 
space vectors; (d) the phase-space vectors are tabulated into 
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non-connected and connected phase-space distribution 
functions; (e) the distribution functions for each sequential 
test state (subset) are compared to the nominal state 
distribution functions via the four measures of dissimilarity 
in Eqs. (7) - (10); (f) each measure of dissimilarity is 
renormalized; (g) these four renormalized measures of 
dissimilarity are used to construct sequential values for the 
composite measure of condition change Q via Eq. (1); and (h) a 

standard least -squares method is used to fit several 
sequential values of C ± to a straight line: 

y s = ai +b (2) 

We choose a window for this straight-line fit with a length of 
m=10 values of C ± (and y ± below) , that is consistent with the 

number of data subsets in each data snapshot. Other values of 
m give inferior indication of condition change. Next, the 

variance, a x 2 # measures the variability of the C s values about 
this straight-line fit: 

<*i 2 = 2, (y, - CjVOn-l) (3) 

Finally, the end-of-life statistical indicator, G, measures 
the variability of next in values of C ± about an extrapolation 
of the straight-line fit as follows: 

a = ( Yi - cjyo* (4) 

The index, i, runs over the m values of C s and y ± in Eqs. (2) - 
(4) . G has the form of a conventional chi-squared statistic, 
but that notation is not used to avoid confusion with the two 
chi-squared PSDMs, U{% 2 ) and Uix c 2 ) • Moreover, this disclosure 
does not use the G measure for a statistical test of null 
hypothesis. Rather in this invention, G is a relative measure, 
providing end-of-life forewarning, plus an indication of the 
failure onset. 

[0014] Other objects and advantages of the invention, besides 
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those discussed above, will be apparent to those of ordinary- 
skill in the art from the description of the preferred 
embodiments, which follows. In the description reference is 
made to the accompanying drawings, which form a part hereof, 
and which illustrate examples of the invention. These 
examples, however are not exhaustive for the various 
embodiments of the invention, and therefore reference is made 
to the claims, which follow the description for determining 
the scope of the invention. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0015] Figs. 1A-1G show graphs of individual three-phase 
currents I ir individual three-phase voltages V it and three- 
phase electrical power vs. time for an electric motor; 

[0016] Figs. 2A-2L show graphs of four phase-space 
dissimilarity measures (rows of subplots) vs. dataset number 
for the three components of three-phase electrical power data 
I S V S from an air-gap seeded- fault in an electrical motor 

(columns of subplots) ; 

[0017] Figs. 3A-3G show graphs of the three components of 
inboard acceleration A it the three components of inboard 
velocity V it and the resultant vibration (mechanical) power; 

[0018] Figs. 4A-4D are graphs of four dissimilarity measures 
vs. the dataset number from the air-gap seeded-fault vibration 
power; 

[0019] Fig. 5 shows the placement of monopolar scalp EEG 
electrodes on the head of a human patient (as viewed by 
looking down on the head from above, with the nose at the top, 
and ears on the left and right sides) , each of which measures 
a local voltage with respect to some reference electrode; 

[0020] Figs. 6A-6C are graphs of three channels of 
acceleration data vs. BIN number for a machine test; 

[0021] Figs. 7A-7L show samples of the three components A s of 
tri -axial acceleration data for the air-gap seeded-fault test 
over successively shorter time intervals, displaying complex, 
nonlinear features; 

[0022] Figs. 8A-8D are graphs of vibration power over 
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successively shorter time intervals, showing complex, 
nonlinear features; 

[0023] Figs. 9A-9F show graphs of phase-space dissimilarity 
measures from vibration power versus time: (a) - (d) the four 
renormalized PSDM;~(e) composite measure, C it of the four PSDM; 

(f ) end-of-life indicator, G (solid) , running maximum of G 

(dashed), and ratio, R, of successive maxima (-.-) in G. 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

[0024] The following description incorporates the methods 
disclosed in U.S. Pat. Appl . No. 10/195,626 filed July 12, 
2002. To the extent that the methods of the present invention 
build upon the methods disclosed there, the disclosure of that 
application is hereby incorporated by reference. 

[0025] The methods of this disclosure can be applied to 
electric motor predictive maintenance, other machinery and 
physical processes, as well as biomedical data for diagnosis 
and treatment of human patients. In one example, data sets 
were recorded in snapshots of 1.5 seconds, sampled at 40 kHz 

(60,000 total time-serial samples), including three-phase 
voltages and currents from an electric motor. The subsequent 
description describes analysis of an induced fault, referred 
to as a seeded fault. 

[0026] In such a motor, the components of the three-phase 
voltages (V s ) and currents (Ij , (see Figs. 1A-1F) can be 

converted into an instantaneous motor power, P = IiV± t where 
the sum runs over the three phases. The electrical power, P, 

(Fig. 1G) is advantageous in forewarning of critical events, 
because it reduces the number of data channels that must be 
analyzed to detect the condition change. The PSDM were 
computed for each of the three-phase currents and voltages 
separately, and compared with the condition change for 
components of the three-phase power. Figs. 2A-2L show 
examples of the four dissimilarity measures for each component 

(IiVj of three-phase power. All four dissimilarity measures 
of I 2 V 2 (middle column) show an almost linear increase, and all 
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four dissimilarity measures of I 3 V 3 increase monotonically. 
All four dissimilarity measures of 1^ increase from test #1 
(no fault) to test #2 (first seeded fault) , then remain 
constant from test #2 (first seeded fault) to test #3 (second 
seeded fault in addition to the first seeded fault) . Thus, it 
was determined that the phase space dissimilarity measures 
(PSDM) method described above can extract condition change 
from one (or more) components of three-phase electrical power 
as well as from individual phase currents and voltages. The 
advantage of selecting electrical power as the process- 
indicative signal is the reduction in the number of channels 
of data for analysis. 

[0027] The data for this same test also included tri-axial 
acceleration data from inboard (IB) and outboard (OB) motor 
locations. Vibration (mechanical) power can be determined from 
tri-axial acceleration, as follows. Acceleration, a, is a 
three-dimensional vector that can be integrated once in time 
to give velocity (vector) , v = J a dt. Mass, m, times 
acceleration (vector) is force (vector) , F = ma. The dot- 
product of force and velocity converts these vector quantities 
into the mechanical power (scalar) , P - F»v = zna»J a dt. 

[0028] Figs. 3A-3G show the three components of inboard 
acceleration (Figs. 3A-3C) , as well as the corresponding three 
components of inboard velocity (Figs. 3D-3F) . The noteworthy 
feature of these plots is the clear periodicity with 
additional complex nonlinear features, which also appear in 
the resulting vibration power (Fig. 3G) . 

[0029] Figs. 4A-4D show a nearly linear rise in all four PSDM 
as computed from the vibration power of the previous 
paragraph. To obtain this result, a sequence of single 
parameter searches was executed: varying S with X fixed, then 
varying X with S set to the value from the first search. This 
approach, yielded S=QS and X=42 for d=3 . Other dimensions 
(d=2 and 4) give inferior results. Subsequently, an exhaustive 
search was made over a range of S and X for d=3 to obtain best 
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condition change at S=84 and Jt=47. 

[0030] The above examples show that electrical power has 
renormalized values of PSDM in the range of 5-55 standard 
deviations from the mean; two of the PSDM rise monotonically 

(not linearly) , thus yielding less desirable change 
discrimination. It was found that vibration power has PSDM 
values on the^ order of unity; all four PSDM increase linearly, 
thus yielding good change discrimination. A search to find 
one (or more) data channels for consistent event forewarning 
requires the analysis of all three components of acceleration, 
or all three components of velocity, or all three components 
of three-phase electrical power, or all six components of 
three-phase voltage and current. However, only one channel of 
electrical power or one of vibration (mechanical) power needs 
to be analyzed to obtain consistent failure forewarning. As 
one normally skilled in the art can readily appreciate, the 
savings in computation effort is an important and novel 
improvement to the art. 

[0031] Another example uses multi-channel EEG data that is 
combined into a multi-channel PS reconstruction for critical- 
event forewarning. Fig. 5 shows an array of monopolar scalp 
EEG electrodes, each of which measures a local voltage with 
respect to some reference electrode. The voltage difference 
between adjacent monopolar electrodes creates a "bipolar 
montage." One bipolar channel is (for example) the difference, 
F4-C4, on the right hemisphere; another is F3-C3 on the left 
hemisphere. Each bipolar channel extracts a time-serial signal 
that potentially allows removal of most locally-generated 
artifacts, while retaining a differential measure of the brain 
activity. Our analysis of EEG data has determined that data 
from the bipolar montage without artifact filtering provides a 
reasonable alternative to artifact filtering, yielding better 
forewarning than the monopolar montage with (or without) 
artifact filtering (and not much different from bipolar 
forewarning with artifact filtering) . Moreover, two (or more) 
of these bipolar channels can be combined via the multi- 
channel phase-space vector, to provide forewarning of 
epileptic seizures that is comparable to the analysis of all 
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nineteen of the individual monopolar channels . As one normally 
skilled in the art can readily appreciate, the sayings in 
computational effort (one two-channel bipolar PS 
reconstruction vs. nineteen individual PS reconstructions) is 
an important and novel improvement. 

[0032] After choosing the appropriate process -indicative data 
in a manner that reduces the analysis effort, as discussed 
above, the present invention next provides an improved method 
to identify inadequate quality data. Raw process -indicative 
data can come . from a variety of sensors, such as 
electroencephalogram data, electrocardiogram data, chest 
sounds, single- or multi-axial accelerometer (s) , single- or 
multi-phase motor current (s) and voltage (s). This data is 
sampled at equal time intervals, x, starting at time t^, 

yielding a sequence of time serial data, Xj = x(t 0 + ix) , i = 
1, 2, N. These signals are digitized and recorded as either 
integer or decimal numbers, which have a range between from 
some minimum value (x^J to maximum value (x^J . 
[0033] Determination of the data quality from the time serial 
data, x lf involves the following steps: (a) sorting the x x 
values into ascending order from x^, to x^; (b) determining 
the number of unique signal values (n) ; (c) obtaining the 
occurrence frequency {F k ) for each unique signal value (v k ) ; 
(d) estimating the occurrence probability for a single outlier 
in a normal distribution by solving: 1/n = # erfc(z/^2) . Here, 
n is the number of unique signal values from step (b) ; erfc is 
the complementary error function; z is the number of standard 
deviations from the mean for such a single outlier; and the 
mean occurrence frequency is estimated as N/n. Next, (e) 
plotting F k versus v k in each BIN, as shown (for example) in 
Figs. 6A-6C; (f) reject any outlier bin with [v k - (N/n)]/o 2 > 
z, where a 3 is the standard deviation of the occurrence 
frequency over all of the unique signal values from step (b) . 
In this example, Figs. 6A-6C show frequency distributions from 
real test data. Huge spikes in log 10 FREQUENCY occur above the 
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underlying baseline (top two subplots for channel 1 (CI) and 
channel 2 (C2), indicating an accumulation of data values in 
one (or a few) unique signal value (s) . These spikes far 
exceed the above limit. Another problem is that the cl- and 
c2-signals have a narrow range of values (-0.099 to +0.089), 
indicating that the signal is improperly scaled for 
digitization. Moreover, the cl- and c2 -distribution functions 
have large variability, which further indicates excessive 
noise. 

[0034] Table 1 summarizes the results for samples of tri- 
axial accelerometer data that was acquired from a test motor. 
The data were sampled for 10 seconds at 128 kHz for a total of 
1,280,000 points from a shaft crack test. The table below 
shows the minimum (min) and maximum (max) values by dataset 
for each of three accelerometer channels (cl - c3) . 

Table 1 . Tri-axial Acceleration Data 

min(cl) max(cl) min(c2) max(c2) min(c3) max(c3) #BINS(c3) 

-.08755 +.07760 -.08748 +.06009 -.00323 -.00254 11 

-.08899 +.07375 -.08247 +.05514 -.00322 -.00261 10 

-.09895 +.07684 -.08893 +.05906 -.00316 -.00261 9 

-.09257 +.08611 -.07705 +.05892 -.00316 0.0 11 

-.08549 +.07725 -.08364 +.05572 -.00309 -.00247 10 

These data have discrete values with five significant figures 
(x.xxxx) for all three channels. The right -most column of 
Table 1 (with the header "#BINS") shows the number of different 
discrete c3 -values, which varies between nine and eleven over 
a very narrow range of signal values (-.0032 to +.0026). This 
small number of discrete signal values is also seen in Fig. 
6C, corresponding to a data precision of lo£r 2 (11) ■ 3.45. This 
low data precision is well below that attainable by modern 
data acquisition systems, which can acquire data with >10 bits 
of precision, and more typically 16 (or more) bits. Thus, the 
c3 -signal is due to bit chatter (noise) , rather than a real 
signal . 

[0035] Rejection of this inadequate quality data from further 
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analysis avoids meaningless or nonsensical results. In the 
prior art represented by U.S. Pat. No. 5,815,413, the data 
range was partitioned into an equal number of non- overlapping 
uniform- size bins to obtain the occurrence frequency. The 
improvement of this disclosure is the use of the unique signal 
values in the raw data. As one of ordinary skill in the art 
can appreciate, this novel improvement uses the explicit 
features in the data to determine its quality, rather than 
artificial constructs of the prior art (such as analyst -chosen 
bins to tabulate the occurrence distribution function) . 
[0036] Data of the types discussed above, may also have 
artifacts that would otherwise confound the analysis, and 
therefore must be removed. In Clapp et al., US Patent No. 
5,626,145, "Method and Apparatus for Extraction of Low- 
Frequency Artifacts from Brain Waves for Alertness Detection," 
a method to remove low- frequency artifacts (such as eye-blinks 
in EEG data or breathing artifacts from ECG data) via a zero- 
phase quadratic filter is disclosed. This filter uses a 
moving window of process -indicative data points, e,, with the 
same number of data points, w, on either side of a central 
point. A quadratic curve is fitted in a least-squares sense to 
these 2w+l data points. The central point of this fit is taken 
as the best estimate of the low-frequency artifact, £,. The 
residual signal, g t - e t - f„ has essentially no low-frequency 
artifact activity. 

[0037] However, data may have both low- frequency artifacts 
(such as eye -blinks in EEG, or breathing in electrocardiogram, 
ECG, data) , and high-frequency artifacts (such as 60-Hz noise 
in ECG) . The novel improvement taught in this disclosure is 
application of the original zero-phase quadratic filter first 
to extract the low- frequency fj-data of the previous paragraph 
using an appropriate (small) value for the filtering -window 
width, w, which effectively removes the high-frequency 
artifacts. Subsequently, the same filter is applied to this f a - 
data with an appropriate (large) filtering window width for w, 
which removes the low- frequency artifacts from the resultant 
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second- stage g^data of the previous paragraph. The use of the 
zero-phase quadratic filter in two stages as applied to PSDM 
data analysis has not been seen in the prior art and improves 
forewarning analysis, for example in ECG data with both low- 
frequency breathing artifacts and high-frequency 60 -Hz noise.. 
As one normally skilled in the art can appreciate, the removal 
of both low- and high-frequency artifacts while preserving the 
nonlinear data features is an important and novel component of 
the analysis. 

[0038] Next, each artifact-filtered value, g i7 is converted 
into a symbolized form, s i9 that is, one of S different 
integers, 0,1, . . . , S-l: 

0 £ Si = INT[S( gi - g min )/(g max - g m±n ) ] < S - 1. (5) 

Here, the function JOT converts a decimal number to the 
closest lower integer; g^ and g TOX denote the minimum and 
maximum values of artifact-filtered data, g it respectively, 
over the baseline (nominal state) data. To maintain S distinct 
symbols, the following expression holds, namely s i = 5.-1 when 
9i = Ska*- Expression (5) creates symbols that are uniformly 
spaced between the minimum and maximum in signal amplitude 
(uniform symbols) . Alternatively, one can use equiprobable 
symbols, by ordering all N base case data points from the 
smallest to the largest value. The first N/S of these ordered 
data values correspond to the first symbol, 0. Ordered data 
values (N/S)+l through 2N/S correspond to the second symbol, 
1, and so on up to the last symbol, 5-1. By definition, 
equiprobable symbols have non-uniform partitions in signal 
amplitude and present the advantage that dynamical structure 
arises only from the phase-space reconstruction. Large 
negative or positive values of g 4 have little effect on 

equiprobable symbolization, but significantly change the 
partitions for uniform symbols. Moreover, the mutual 
information function is a smooth function of the 
reconstruction parameters for equiprobable symbols, but is a 
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noisy function of these same parameters for uniform symbols. 
Thus, equiprobable symbols provide better indication of 
condition change than uniform symbols when constructing a 
(non) connected phase space . 

[0039] Phase-space (PS) construction uses time-delay vectors, 

= [s iy s iA / • • . / Si+uf.i)X] to unfold the underlying 

dynamics. Critical parameters in this approach are the time 

delay, X, and system dimensionality, d, and the type and 

number of symbols, S. However, information exchange between 
(for example) complex machinery components or between 
different brain locations connect the local processes. This 
disclosure provides a novel method to combine several channels 
of time-serial data into a multi-channel phase-space (PS) 
reconstruction, which has the form of a dC-dimensional vector: 

Sj=[s(l) it S(l) i+ X ,~ / 8(l)i+(d-l)Xf- / s(Oif 8(C) i4 X,~ s(C)i+(d-i)d • 

Here, Sj(j) denotes the symbolized data for j-th channel at 
time t £ . C denotes the total number of data channels in the 
multi-channel phase-space reconstruction. 

[0040] Conversion of the multi-channel time serial data to a 
sequence of phase-space vectors via the method of the previous 
paragraph divides the phase space into SF* discrete bins, which 

are uniquely labeled via base-S arithmetic: 



The resulting distribution function (DF) is a discretized 
density, which is obtained by counting the number of points 
that occur in each phase space bin. The population of the ith 
bin of the distribution function, is denoted Qx, for the base 
case, and R± for a test case, respectively. 

[0041] The test case is compared to the base case via 
dissimilarity measures between Qj with R £ that are defined by: 



c 



d-1 



0 <, yU) = yist) = S S By),* S^- 1 ' £ S" - 1 
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Here, the summations run over all of the populated phase space 
cells. These measures account for changes in the geometry and 
visitation frequency of the resultant geometric object, called 
an "attractor. 1 1 The % 2 measure is one of the most powerful, 
robust, and widely used statistics for comparison between 
observed and expected frequencies. In this context, % 2 is a 
relative measure of dissimilarity, rather than an unbiased 
statistic for accepting or rejecting a null statistical 
hypothesis. The L distance is the natural metric for 
distribution functions by its direct relation to the total 
invariant measure on the phase-space attractor. Consistent 
calculation of these measures obviously requires the same 
number of points in both the base case . and the test case 
distribution functions, identically sampled; otherwise, the 
distribution functions must be properly rescaled. 
[0042] The process dynamics yield successive PS points, as 
prescribed by -» Thus, a discrete representation of 

the process flow is obtained in the form of a 2 do dimensional 
multi- channel PS vector, Sj = [s^, s^J , that is formed by 
adjoining two successive vectors from the dC- dimensional 
reconstructed PS. Sj is a vector for the multi-channel 
connected phase space (CPS) . As before, Q and R denote the CPS 

distribution functions for the base case and test case, 
respectively. The measure of dissimilarity between the two 
distribution functions for the CPS, signified by the w c" 
subscript are thus defined as follows: 



The first subscript in Eqs. (9) -(10) denotes the initial PS 
point, and the second subscript denotes the sequel PS point in 
the CPS pair.. These CPS measures have higher discriminating 
power than unconnected PS measures of dissimilarity. Indeed 




(9) 



4 = Sa-*l- 



(10) 
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the measures defined in Eqs. (7) -(10) satisfy the following 
inequalities: J? <> L, Xc * L o L * L o and f < X* * where Xc and 
L c are dissimilarity measures for connected phase space and 
and L are dissimilarity measures for non-connected PS. 
[0043] To assure robustness, the construction of the base 
case data requires careful statistics to eliminate outlier 
base case subsets. For a description of this part of the 
methodology, reference is made to Hively, U.S. Pat. Appl . No. 
10/195,626 filed July 12, 2002. 

[0044] A consistent method of comparison is needed to 
interpret the disparate range and variability of the 
dissimilarity measures. To this end, the dissimilarity 
measures are renormalized, as described next. The B non- 
outlier base case subsets from the previous paragraph are 
compared to each test case subset, to obtain the corresponding 
average dissimilarity value, V if of the ith subset for each 
dissimilarity measure. Here, V denotes each dissimilarity 
measure from the set, V = {L, L c , jf, and Xc) • The mean value, 
V, and the corresponding standard deviation, a 2 , of the 
dissimilarity measure V are calculated by comparison of all of 
the unique pairs of base case subsets, after the outliers have 
been eliminated. The renormalized dissimilarity is the number 
of standard deviations that the test case deviates from the 
base case mean: U(V) m \v ± - v\/c 2 . 

[0045] Once the renormalized measures for the test and base 
cases have been obtained, a threshold, U c , is selected for each 
renormalized measure U to distinguish between normal (base) 
and possibly abnormal (test) regimes. The choice of a 
reasonable threshold is critical for obtaining robust, 
accurate, and timely results. A forewarning indication is 
obtained when a renormalized measure of dissimilarity exceeds 
the threshold, U £ U c , for a specified number, N^, of 
sequential occurrences. 

[0046] The remainder of the PSDM analysis is described by 
Hively, U.S. Pat. Appl. No. 10/195,626 filed July 12, 2(M)2, 
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cited above. As one normally skilled in the art can readily 
appreciate, the novel improvement to the prior art involves 
combining multiple -channel time-serial data into a one multi- 
channel time-delay phase-space (PS) reconstruction, rather 
than performing the PS reconstruction for each of several 
channels . 

[0047] An alternative embodiment of this invention is 
illustrated by the following example. The vibration power 
data, as discussed above, was acquired in 10- second snapshots 
at 50 kHz at fifteen-minute intervals; see Table 2 for the 
interval (Ax) between each snapshot. Each 500,000-point 
snapshot was divided into ten non-overlapping subsets of 
50, 000-points each for subsequent analysis. This data 
characterized an accelerated- failure test of an over-loaded 
gearbox. The test protocol involves a break- in period at the 
nominal (IX) load (530 ft -lbs) for one hour, followed by twice 
(2X) or three times (3X) the normal load, as shown in Table 2 
The end-of-life failure included pinion gear damage, broken 
teeth, and a sheared shaft . The failure occurred after 162.5 
hours, corresponding to 650 snapshots. Figs. 7A-7L show the 
very complex, nonlinear features in the three components of 
the tri- axial accelerometer data (a) . The individual snapshots 
were combined into one long dataset (12.7 GB) . Acceleration 
was subsequently converted to a long stream (4.1 GB) of 
vibration power (Figs. 8A-8D) , using the previously described 
method to convert tri -axial acceleration into power via time- 
integration to velocity (v = Jadt) with a subsequent vector 
dot -product to produce power (P ~ a» v) . Figs. 9A-9D show the 
PSDM from this data. The phase-space parameters are N=50,000, 
5=274, d=2, JL»1. All four dissimilarity measures rise 
systematically (Figs. 9A - 9D) to provide forewarning of the 
failure. However, a more robust and specific end-of-life (EOL) 
indicator is needed. We observe that all four of the PSDMs 
have similar tends, suggesting the definition of a composite 
measure, C 1# as the sum of the four renormalized PSDMs for the 
i-th dataset (Fig. 9E) : 
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C, = Utf) + tf( Xc 2 ) + U(L) + E7(L C ) (11) 

This composite measure is expected to be more robust than any 
one of the individual measures of phase space dissimilarity, 
while accurately indicating condition change. Subsequent 
description presents another novel improvement of this 
disclosure to provide end-of-life forewarning, and indication 
of failure onset. 

[0048] The end-of-life indication from the composite measure 
of the previous paragraph is quantified as follows. We use 
contiguous, non-overlapping windows of C ± to obtain the best 
straight-line fit in a least-squares sense: 

Yi = ai +b (12) 

The window length of m=10 values of C ± (and y t below) is chosen 
consistent with the number of subsets in each snapshot. Other 
values of m give inferior indication of condition change. 

Next, the variance, a x 2 , measures the variability of the C ± 
values about this straight-line fit: 

a x 2 = 2, (y, - cyVUn-U (13) 

[0049] Other fits (such as, quadratic, cubic, and quartic) 
give inferior condition change indication, due to the poor 
extrapolation outside the fitting window. Finally, G measures 
the variability of next m values of C, about an extrapolation 
of this straight -line: 

G = Z 4 (K - CjVa, 2 (14) 

The index, i, in Eqs. (12) - (14) runs over the in values of C x 
and y x . Note that G has the form of a conventional chi-squared 
statistic, but we do not use that notation to avoid confusion 
with the two % 2 PSDMs, C7(x 2 ) and £7(x c 2 ) - A statistical test for 
G would involve (for example) the null hypothesis that 
deviations from the straight-line fit are normally 
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distributed. Standard chi-squared statistical tables give the 
corresponding value of G < 28.5 for m = 10 degrees of freedom 
with a probability of one out of the number (650) of snapshots 
(1/650 ~ 1.5 x 10" 3 ) . However, we observe many instances with 
G>28.5 throughout this test sequence. These outliers occur 
because the underlying tri- axial acceleration has dynamical 
correlations, thus violating the requirement for independent, 
identically distributed samples. Instead, we use G (solid 
curve in Fig. 9F) as a relative measure of end- of -life. As 
time moves forward, we obtain the running maximum of G, as the 
larger of G and (dashed curve in Fig. 9F) . We neglect the 

first six G-values (for example) to avoid startup transients. 
This running maximum rises in modest increments to 376 over 
the first 159.75 hours of the test, while intermediate values 
of G fall well below G^. The chain curve (-.-) in Fig. 9F is 
the ratio, R = (G^J J (G^) k _ 2 , of the current running maximum 
in G, (GfcaJjt, to the previous maximum in G, (Gj^. G m rises 
to 2,493 at 160 hours, with a corresponding R = 6.62, which is 
well above the largest non-end-of-lif e value, R = 2.22, at 2 
hours. Thus, G provides clear end-of-life forewarning, plus 
indication of the failure (G = 244,655). Table 2 gives: (a) 
the largest non-EOL value of R (iWJ and the corresponding 
value of G (G m ) ; (b) values of R {R m ) and G (G TOL ) that 
indicate the end of life, and the matching time (T^/T^) ; (c) 
the value of G at failure onset (G^^) and the corresponding 
time (Tann/T^n); and (d) the f ailure-endpoint time (T FAIL ) . 



Table 2 • Summary of Gearbox Overload Results 
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Over- 
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At 

min. 
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RxOL 


G E0L 
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hr. 


36 


2X 


15 


2.22 


376 


6.62 


2,493 


.985 


244,655 


0.998 


162.50 


37 


3X 


l 


1.79 


333 


B.07 


2,690 


.956 


16,284 


0.996 


8.55 


38 


3X 


l 
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11.71 


13,486 


.938 


48,379 


0.990 


4.02 


39 


2X 


i 


2.32 


B53 


3.89 


5,231 


.980 


5,231 


0.980 


8.60 
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1" 


3X 


1 


2.88 


1,151 


29-03 


33,415 


.972 


44,552 


0.994 


8.60 | 



Runs #36-38 have largest non-EOL values: R^o h = 6.20 and G NEOL « 
376. The smallest EOL values are: R^ = 6.62 and = 2,493. 

Thus, limits (for example) of > 6.4 and G > 1,800 provide 
EOL forewarning. Moreover, we find that the largest EOL value 
of G m = 13,486, while the smallest failure-onset value is 
g onsbt = 16,284. Thus, an intermediate value (for example) of G 
> 15,000 distinguishes the EOL from failure onset forewarning. 
This approach gives quantitative limits for transitions from 
nominal operation (green-light in a traffic signal metaphor) , 
to forewarning of failure (yellow light for caution) , and 
finally to failure onset (red-light for stop) . 

[0050] Run #39 involves a different test protocol: a one-hour 
break- in period at nominal load (IX) , followed by 2X load for 
two hours, after which the load alternates between 3X and 2X 
loads for ten and five minutes, respectively. Run #39 seeks 
failure forewarning in the presence of load changes. Table 2 
shows that the above limits for G and R also distinguish 
between the non-EOL (green) and EOL (yellow) states for the 
3X-portion of this test, because the higher overload drives 
the failure. These limits do not apply to the 2X test, due to 
the reduced damage at the lower overload. Unsurprisingly, a 
different limit of G > 38,000 (for example) distinguishes 
between the EOL and failure onset f orewarnings , due to the 
change in test protocol. The green-yellow-red traffic-signal 
metaphor still applies for this test. 

[0051] As one of ordinary skill in the art can appreciate, 
these examples give quantitative limits for transitions from 
normal operation (normal condition) , to forewarning of failure 

(caution condition) , and finally to failure onset (critical 
event), as a novel improvement in the art. 
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[0052] This has been a description of detailed examples of 
the invention. These examples illustrate the technical 
improvements, as taught in the present invention, it will be 
apparent to those of ordinary skill in the art that certain 
modifications might be made without departing from the scope 
of the invention, which is defined by the following claims. 
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CLAIMS 

I claim: 



1. A method for processing data to provide a forewarning 
of a critical event, the method comprising: 

acquiring a plurality of sets of data for at least one 
channel of data for at least one test subject or process; 

computing a renormalized measure of dissimilarity from 
distribution functions derived from a phase space for each 
respective channel of data; 

comparing said renormalized measure of dissimilarity to a 
threshold (C7 C ) for a number of occurrences {Nocc) to indicate a 
condition change in said renormalized measure of 
dissimilarity; 

and further characterized by: 

detecting a simultaneous condition change in a plurality 
(N S im) of renormalized measures of dissimilarity to determine a 
forewarning of the critical event; and 

wherein said one channel of data corresponds to a 
parameter that is calculated from a plurality of parameters 
corresponding to a plurality of channels of data. 

2. The method of claim 1, wherein the test subject is a 
human patient. 

3. The method of claim 1, wherein the test subject is a 
mechanical device or physical process. 

4. The method of claim 1, wherein the process -indicative 
data is three-phase electrical power. 

5. The method of claim 1, wherein the process- 
indicative data is vibration mechanical power. 

6. The method of claim 1, wherein the process -indicative 
data is a difference between two channels of EEG data. 
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7. The method of claim 1, further characterized by: 
performing a first filtering of each set of data with a 

zero-phase quadratic filter that filters out high-frequency 
artifacts; and 

performing a second filtering of each set of data with a 
zero-phase quadratic filter to filter out low-frequency 
artifacts. 

8. The method of claim 1, further characterized by: 
sorting the data values into ascending order from a 

minimum to a maximum; 

determining the number of unique signal values (n) and 

the corresponding relative occurrence frequency (F k ) for each 

unique signal value (v k ) ; 

displaying a graph of frequency (F k ) versus values (v k ) in 
each bin in a connected phase space; and 

discarding data that has [v k - (N/n) ] /a 3 > z, where the 
value of z is determined by solving l/n = M erfc(z/V2) , and a 3 
is the standard deviation in the occurrence frequency. 

9. The method of claim 1, with an alternative embodiment 
for event forewarning, characterized by determining a sequence 
of renormalized phase space dissimilarity measures from data 
sets for the test subject or process; summing said 
renormalized measures into a composite measure, C i# for the i- 
th data set; performing a least-squares analysis over a window 
of m points of the said composite measure to obtain a straight 
line, y^ai+Jb, that best fits said composite data in a least- 
squares sense; determining the variance, a x 2 = E ± (y* - d) 2 /(m- 
1) , of said composite measure with respect to the straight 
line fit; obtaining the variability of the sequel window of m 

sequential points via the statistic, G = 2± (y± - Ci) 2 /<Ti 2 ; 
comparing said value of G to the running maximum value of the 
same statistic, G^*; determining the forewarning of or failure 
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onset of a critical event (such as a machine failure) , when G 
is significantly more than Gmax; obtaining the ratio, R = 
(GmaJk/ (3nax)*-i' of the present and previous running maximum in 
G; and determining the forewarning of a critical event when R 
is significantly more than some limit. 

10. A method for processing data to provide a 
forewarning of a critical event, comprising: 

acquiring a plurality of sets of data for at least two 
channels of data for at least one test subject or process; 
and characterized by: 

computing to a multi-channel time-delay phase-space (PS) 
construction, which has the form: y(i) = [s(l) £ , s(l) i+ x , 
s(1) U2 a,- , 3(2) , s(2) t ^, s(2) x , 2 x, ~, s(c) if s(c)^ f s(c) M x, 
where s(c) denotes the symbolized data for c-th channel; 

computing a renormalized measure of dissimilarity from 
distribution functions derived from a non-connected phase 
space for the multi-channel of data; 

comparing said renormalized measure of dissimilarity to a 
threshold (C7 C ) for a number of occurrences {Nocc) to indicate a 
condition change in said renormalized measure of 
dissimilarity; and 

detecting a simultaneous condition change in a plurality 
(N S im) of renormalized measures of dissimilarity to determine a 
forewarning of the critical event. 

11. The method Of claim 10, further characterized by: 
performing a first filtering of each set of data with a 

zero-phase quadratic filter that filters out high-frequency 
artifacts; and 

performing a second filtering of each set of data with a 
zero-phase quadratic filter to filter out low- frequency 
artifacts . 

12. The method of claim 10, using an alternative 
embodiment for event forewarning, further characterized by 
determining a sequence of renormalized phase space 
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dissimilarity measures from data sets collected during 
increasingly severe fault conditions; summing said 
renormalized measures into a composite measure, C ir for the i- 
th data set; performing a least -squares analysis over a window 
of m points of the said composite measure to obtain a straight 
line, y—ai+b, that best fits said composite data in a least- 
squares sense; determining the variance, a x 2 = S± (yi - Ci) 2 /{m- 
1), of said composite measure with respect to the straight 
line fit; obtaining the variability of a sequel window of m 

sequential points via the statistic, G = Z ± (y± - C±) 2 /oi 2 ; 
, comparing said value of 6 to the running maximum value of the 
same statistic, Gma X ; and determining the onset of a critical 
event, such as forewarning of a machine failure, when G is 
significantly more than G (non-end-of -lif e) or when R is 

significantly more than £(non-end-of -lif e) , or detection of 
failure onset when G is significantly greater than G(end-of- 
life) . 

13. The method of claim 10, wherein the test subject is 
a human patient. 

14. The method of claim 10, wherein the test subject is 
a mechanical device or physical process. 

15. The method of claim 10, further characterized by: 
sorting the data values into ascending order from a 

minimum to a maximum; 

determining the number of unique signal values (n) and 

the corresponding relative occurrence frequency {F k ) for each 

unique signal value (v k ) ; 

displaying a graph of frequency {F k ) versus values (v*) in 
each bin in a connected phase space; and 

discarding data that has [v k - (N/n) ] /a 3 > z, where the 

value of z is determined by solving l/n = & erfc(z/V2), and o 2 
is the standard deviation in the occurrence frequency. 
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16. The method of claim 10, wherein the process- 
indicative data is three-phase electrical power. 

17. The method of claim 1, wherein the process- 
indicative data is vibration mechanical power. 

18. The method of claim 1, wherein the process- 
indicative data is a difference between two channels of EEG 
data. 
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